week 2: linear model and causal inference
geocentric models
Annoucements and such
- We will start using
brmstoday! Install this package now, if you haven’t already.
install.packages(c("brms","tidybayes"))- No office hours Friday
- Lori Olsen
- Projects
Workspace setup
library(here)
library(tidyverse)
library(cowplot)
library(brms)
library(tidybayes)
library(patchwork)# here's a bunch of stuff I include at the beginning of each qmd file.
# it's mostly setting the color palette for figures, but it also controls the default code chunk option (showing you, i.e., echo) and setting the figure quality to ensure that it looks nice
knitr::opts_chunk$set(fig.retina=3, echo=TRUE)
theme_set(theme_cowplot())
default_palettes <- list(
c("#5e8485" , "#0f393a") ,
c("#1c5253" , "#5e8485" , "#0f393a") ,
# palette with 5 colours
c( "#1c5253" , "#e07a5f", "#f2cc8f" , "#81b29a" , "#3d405b" ) ,
# same palette interpolated to 8 colours
c( "#1c5253" , "#e07a5f", "#f2cc8f" , "#81b29a" , "#3d405b" , "#a7a844" , "#69306d" )
)
options(ggplot2.discrete.fill = default_palettes,
ggplot2.discrete.colour = default_palettes)
Workflow
State a clear question.
Sketch your causal assumptions.
Use the sketch to define a generative model.
Use the generative model to build an estimator.
Profit.
Model “recipes”
- Recognize a set of variables to work with. (Data and parameters.)
- Define each variable either in terms of the other variables OR in terms of a probability distribution.
- The combination of variables and their probability distributions defines a joint generative model that can be used to simulate hypothetical observations and analyze real ones.
Here’s an example:
\[\begin{align*} y_i &\sim \text{Normal}(\mu_i,\sigma) \\ \mu_i &= \beta x_i \\ \beta &\sim \text{Normal}(0,10) \\ \sigma &\sim \text{Exponential}(1) \\ x_i &\sim \text{Normal}(0,1) \\ \end{align*}\]
Model for globe-tossing
Here’s the model for last week’s globe-tossing experiment:
\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]
- \(W\) is the observed count of water.
- \(N\) is the total number of tosses.
- \(p\) is the proportion of water on the globe.
The whole model can be read as:
The count \(W\) is distributed binomially with sample size \(N\) and probability \(p\). The prior for \(p\) is assumed to be uniform between 0 and 1.
Model for globe-tossing
Here’s the model for last week’s globe-tossing experiment:
\[\begin{align*} W &\sim \text{Binomial}(N,p) \\ p &\sim \text{Uniform}(0,1) \\ \end{align*}\]
Estimating the posterior using brms
Last week, we used grid approximation to estimate the posterior distribution. In the video you watched for today, McElreath moves on to something called QUADRATIC APPROXIMATION. It’s good to understand what that’s doing, but you and I are moving right along to MCMC. We won’t go into the details of what’s happening for a few weeks, but let’s start with the code to estimate the model.
Let’s say we tossed the globe 9 times and observed 6 waters:
m1 <-
brm(data = list(w = 6), # <1>
family = binomial(link = "identity"), # <2>
w | trials(9) ~ 0 + Intercept, # <3>
prior(uniform(0, 1), class = b), # <4>
iter = 5000, warmup = 1000, seed = 3, chains=1, # <5>
file = here("files/models/m21.1")) # <6>- Data can be a data frame or a list. Just make sure the variable names match what’s in your formula.
- How you assume your outcome variable is distributed.
- The formula for your outcome.
- Priors for every parameter in your model. Here, we only have one parameter, so we only need 1 prior. This happens to be a flat prior between 0 and 1.
- Some choices about how we want our model to run. We’ll go more into this later.
- These models can take a long time to run. You can automatically save the output to a file; when you do this, the next time you run this code, it won’t actually estimate the model, but will instead pull the output from your stated file. Be WARNED: if you change the data or the model code, it will NOT restimate your model until you delete the file.
Sampling from the posterior
Grid approximation gave us the calculated probability of each possible value of our parameter, \(p\). But our method of conducting bayes will no longer give us such a neat solution. Here’s how you get the posterior distribution for \(p\):
samples_from_post = as_draws_df(m1)
samples_from_post## # A draws_df: 4000 iterations, 1 chains, and 3 variables
## b_Intercept lprior lp__
## 1 0.46 0 -2.1
## 2 0.49 0 -1.9
## 3 0.56 0 -1.5
## 4 0.66 0 -1.3
## 5 0.71 0 -1.3
## 6 0.62 0 -1.3
## 7 0.68 0 -1.3
## 8 0.55 0 -1.5
## 9 0.63 0 -1.3
## 10 0.67 0 -1.3
## # ... with 3990 more draws
## # ... hidden reserved variables {'.chain', '.iteration', '.draw'}samples_from_post %>%
ggplot(aes(x=b_Intercept)) +
geom_density(fill = "grey", color = "white") +
labs(x="Proportion water")
Posterior predictive distribution
ppd = posterior_predict(m1)
dim(ppd)## [1] 4000 1ppd## [,1]
## [1,] 6
## [2,] 4
## [3,] 6
## [4,] 5
## [5,] 6
## [6,] 6
## [7,] 6
## [8,] 4
## [9,] 7
## [10,] 7
## [11,] 5
## [12,] 7
## [13,] 8
## [14,] 7
## [15,] 8
## [16,] 9
## [17,] 7
## [18,] 9
## [19,] 8
## [20,] 9
## [21,] 7
## [22,] 9
## [23,] 8
## [24,] 7
## [25,] 6
## [26,] 6
## [27,] 3
## [28,] 4
## [29,] 5
## [30,] 3
## [31,] 5
## [32,] 8
## [33,] 4
## [34,] 7
## [35,] 3
## [36,] 2
## [37,] 5
## [38,] 5
## [39,] 5
## [40,] 5
## [41,] 6
## [42,] 8
## [43,] 7
## [44,] 6
## [45,] 6
## [46,] 6
## [47,] 7
## [48,] 6
## [49,] 7
## [50,] 9
## [51,] 7
## [52,] 3
## [53,] 6
## [54,] 7
## [55,] 4
## [56,] 7
## [57,] 5
## [58,] 2
## [59,] 6
## [60,] 6
## [61,] 7
## [62,] 6
## [63,] 4
## [64,] 6
## [65,] 6
## [66,] 6
## [67,] 8
## [68,] 7
## [69,] 6
## [70,] 4
## [71,] 6
## [72,] 5
## [73,] 8
## [74,] 5
## [75,] 5
## [76,] 3
## [77,] 2
## [78,] 6
## [79,] 9
## [80,] 5
## [81,] 6
## [82,] 7
## [83,] 5
## [84,] 5
## [85,] 6
## [86,] 7
## [87,] 5
## [88,] 5
## [89,] 6
## [90,] 6
## [91,] 6
## [92,] 3
## [93,] 5
## [94,] 5
## [95,] 5
## [96,] 6
## [97,] 7
## [98,] 7
## [99,] 7
## [100,] 8
## [101,] 5
## [102,] 7
## [103,] 8
## [104,] 6
## [105,] 6
## [106,] 8
## [107,] 5
## [108,] 6
## [109,] 6
## [110,] 3
## [111,] 5
## [112,] 4
## [113,] 5
## [114,] 8
## [115,] 9
## [116,] 3
## [117,] 6
## [118,] 5
## [119,] 7
## [120,] 7
## [121,] 8
## [122,] 3
## [123,] 5
## [124,] 7
## [125,] 7
## [126,] 7
## [127,] 5
## [128,] 6
## [129,] 5
## [130,] 6
## [131,] 9
## [132,] 8
## [133,] 4
## [134,] 5
## [135,] 4
## [136,] 8
## [137,] 4
## [138,] 7
## [139,] 7
## [140,] 5
## [141,] 9
## [142,] 5
## [143,] 6
## [144,] 5
## [145,] 6
## [146,] 8
## [147,] 6
## [148,] 6
## [149,] 5
## [150,] 7
## [151,] 6
## [152,] 9
## [153,] 6
## [154,] 7
## [155,] 8
## [156,] 4
## [157,] 6
## [158,] 7
## [159,] 5
## [160,] 3
## [161,] 6
## [162,] 3
## [163,] 5
## [164,] 8
## [165,] 6
## [166,] 6
## [167,] 3
## [168,] 6
## [169,] 7
## [170,] 9
## [171,] 5
## [172,] 5
## [173,] 5
## [174,] 8
## [175,] 7
## [176,] 3
## [177,] 9
## [178,] 9
## [179,] 6
## [180,] 7
## [181,] 5
## [182,] 7
## [183,] 9
## [184,] 7
## [185,] 6
## [186,] 5
## [187,] 8
## [188,] 6
## [189,] 6
## [190,] 5
## [191,] 8
## [192,] 8
## [193,] 6
## [194,] 7
## [195,] 4
## [196,] 5
## [197,] 6
## [198,] 7
## [199,] 8
## [200,] 8
## [201,] 9
## [202,] 4
## [203,] 8
## [204,] 7
## [205,] 3
## [206,] 6
## [207,] 4
## [208,] 6
## [209,] 7
## [210,] 7
## [211,] 7
## [212,] 4
## [213,] 9
## [214,] 5
## [215,] 8
## [216,] 8
## [217,] 9
## [218,] 6
## [219,] 7
## [220,] 5
## [221,] 6
## [222,] 4
## [223,] 3
## [224,] 6
## [225,] 4
## [226,] 3
## [227,] 4
## [228,] 9
## [229,] 6
## [230,] 6
## [231,] 8
## [232,] 7
## [233,] 5
## [234,] 5
## [235,] 7
## [236,] 6
## [237,] 9
## [238,] 7
## [239,] 5
## [240,] 3
## [241,] 3
## [242,] 6
## [243,] 7
## [244,] 8
## [245,] 6
## [246,] 6
## [247,] 1
## [248,] 8
## [249,] 4
## [250,] 5
## [251,] 6
## [252,] 8
## [253,] 7
## [254,] 9
## [255,] 8
## [256,] 7
## [257,] 4
## [258,] 7
## [259,] 7
## [260,] 4
## [261,] 7
## [262,] 8
## [263,] 5
## [264,] 8
## [265,] 4
## [266,] 8
## [267,] 5
## [268,] 7
## [269,] 6
## [270,] 8
## [271,] 6
## [272,] 6
## [273,] 9
## [274,] 5
## [275,] 5
## [276,] 4
## [277,] 9
## [278,] 8
## [279,] 7
## [280,] 7
## [281,] 6
## [282,] 3
## [283,] 5
## [284,] 1
## [285,] 2
## [286,] 1
## [287,] 1
## [288,] 9
## [289,] 6
## [290,] 8
## [291,] 7
## [292,] 9
## [293,] 9
## [294,] 7
## [295,] 4
## [296,] 6
## [297,] 8
## [298,] 6
## [299,] 1
## [300,] 2
## [301,] 5
## [302,] 7
## [303,] 6
## [304,] 7
## [305,] 7
## [306,] 6
## [307,] 8
## [308,] 7
## [309,] 6
## [310,] 5
## [311,] 4
## [312,] 7
## [313,] 7
## [314,] 4
## [315,] 6
## [316,] 8
## [317,] 5
## [318,] 7
## [319,] 5
## [320,] 4
## [321,] 4
## [322,] 4
## [323,] 4
## [324,] 3
## [325,] 5
## [326,] 7
## [327,] 5
## [328,] 5
## [329,] 4
## [330,] 8
## [331,] 4
## [332,] 7
## [333,] 8
## [334,] 5
## [335,] 8
## [336,] 8
## [337,] 8
## [338,] 5
## [339,] 2
## [340,] 3
## [341,] 5
## [342,] 6
## [343,] 7
## [344,] 7
## [345,] 7
## [346,] 4
## [347,] 8
## [348,] 9
## [349,] 8
## [350,] 7
## [351,] 8
## [352,] 6
## [353,] 4
## [354,] 4
## [355,] 8
## [356,] 7
## [357,] 3
## [358,] 1
## [359,] 3
## [360,] 6
## [361,] 4
## [362,] 4
## [363,] 7
## [364,] 4
## [365,] 6
## [366,] 4
## [367,] 7
## [368,] 4
## [369,] 9
## [370,] 8
## [371,] 6
## [372,] 3
## [373,] 3
## [374,] 3
## [375,] 2
## [376,] 3
## [377,] 4
## [378,] 8
## [379,] 2
## [380,] 6
## [381,] 9
## [382,] 6
## [383,] 9
## [384,] 7
## [385,] 7
## [386,] 7
## [387,] 6
## [388,] 5
## [389,] 3
## [390,] 7
## [391,] 7
## [392,] 7
## [393,] 7
## [394,] 4
## [395,] 8
## [396,] 9
## [397,] 5
## [398,] 5
## [399,] 8
## [400,] 7
## [401,] 6
## [402,] 6
## [403,] 7
## [404,] 8
## [405,] 4
## [406,] 6
## [407,] 6
## [408,] 7
## [409,] 9
## [410,] 8
## [411,] 8
## [412,] 7
## [413,] 6
## [414,] 4
## [415,] 7
## [416,] 6
## [417,] 4
## [418,] 4
## [419,] 7
## [420,] 4
## [421,] 5
## [422,] 7
## [423,] 2
## [424,] 6
## [425,] 8
## [426,] 6
## [427,] 3
## [428,] 8
## [429,] 6
## [430,] 5
## [431,] 6
## [432,] 3
## [433,] 4
## [434,] 7
## [435,] 7
## [436,] 8
## [437,] 7
## [438,] 3
## [439,] 3
## [440,] 8
## [441,] 6
## [442,] 8
## [443,] 7
## [444,] 7
## [445,] 7
## [446,] 4
## [447,] 5
## [448,] 5
## [449,] 4
## [450,] 7
## [451,] 3
## [452,] 5
## [453,] 6
## [454,] 2
## [455,] 6
## [456,] 3
## [457,] 8
## [458,] 7
## [459,] 9
## [460,] 6
## [461,] 8
## [462,] 7
## [463,] 8
## [464,] 5
## [465,] 8
## [466,] 7
## [467,] 5
## [468,] 4
## [469,] 6
## [470,] 5
## [471,] 5
## [472,] 8
## [473,] 6
## [474,] 5
## [475,] 6
## [476,] 7
## [477,] 4
## [478,] 5
## [479,] 7
## [480,] 6
## [481,] 7
## [482,] 5
## [483,] 6
## [484,] 6
## [485,] 5
## [486,] 6
## [487,] 5
## [488,] 4
## [489,] 6
## [490,] 6
## [491,] 8
## [492,] 4
## [493,] 2
## [494,] 3
## [495,] 5
## [496,] 7
## [497,] 6
## [498,] 5
## [499,] 3
## [500,] 3
## [501,] 4
## [502,] 5
## [503,] 7
## [504,] 5
## [505,] 4
## [506,] 6
## [507,] 8
## [508,] 6
## [509,] 6
## [510,] 7
## [511,] 8
## [512,] 6
## [513,] 6
## [514,] 6
## [515,] 6
## [516,] 6
## [517,] 6
## [518,] 7
## [519,] 4
## [520,] 4
## [521,] 2
## [522,] 1
## [523,] 3
## [524,] 5
## [525,] 4
## [526,] 8
## [527,] 8
## [528,] 5
## [529,] 4
## [530,] 7
## [531,] 7
## [532,] 9
## [533,] 8
## [534,] 5
## [535,] 6
## [536,] 5
## [537,] 6
## [538,] 7
## [539,] 6
## [540,] 7
## [541,] 3
## [542,] 5
## [543,] 6
## [544,] 7
## [545,] 7
## [546,] 8
## [547,] 5
## [548,] 6
## [549,] 6
## [550,] 4
## [551,] 8
## [552,] 7
## [553,] 3
## [554,] 6
## [555,] 6
## [556,] 4
## [557,] 3
## [558,] 3
## [559,] 3
## [560,] 4
## [561,] 4
## [562,] 6
## [563,] 6
## [564,] 5
## [565,] 6
## [566,] 4
## [567,] 5
## [568,] 7
## [569,] 8
## [570,] 7
## [571,] 7
## [572,] 2
## [573,] 4
## [574,] 4
## [575,] 5
## [576,] 8
## [577,] 6
## [578,] 7
## [579,] 5
## [580,] 6
## [581,] 5
## [582,] 5
## [583,] 6
## [584,] 6
## [585,] 6
## [586,] 5
## [587,] 2
## [588,] 5
## [589,] 5
## [590,] 2
## [591,] 5
## [592,] 5
## [593,] 7
## [594,] 6
## [595,] 4
## [596,] 6
## [597,] 6
## [598,] 4
## [599,] 5
## [600,] 8
## [601,] 7
## [602,] 4
## [603,] 5
## [604,] 7
## [605,] 6
## [606,] 8
## [607,] 8
## [608,] 2
## [609,] 8
## [610,] 6
## [611,] 2
## [612,] 9
## [613,] 5
## [614,] 5
## [615,] 5
## [616,] 3
## [617,] 4
## [618,] 8
## [619,] 8
## [620,] 7
## [621,] 5
## [622,] 7
## [623,] 6
## [624,] 7
## [625,] 3
## [626,] 4
## [627,] 6
## [628,] 9
## [629,] 4
## [630,] 4
## [631,] 3
## [632,] 7
## [633,] 8
## [634,] 5
## [635,] 5
## [636,] 5
## [637,] 6
## [638,] 0
## [639,] 0
## [640,] 3
## [641,] 5
## [642,] 5
## [643,] 6
## [644,] 7
## [645,] 5
## [646,] 5
## [647,] 3
## [648,] 4
## [649,] 7
## [650,] 5
## [651,] 8
## [652,] 6
## [653,] 7
## [654,] 2
## [655,] 4
## [656,] 7
## [657,] 9
## [658,] 8
## [659,] 6
## [660,] 8
## [661,] 7
## [662,] 5
## [663,] 4
## [664,] 5
## [665,] 8
## [666,] 7
## [667,] 9
## [668,] 8
## [669,] 7
## [670,] 5
## [671,] 7
## [672,] 7
## [673,] 4
## [674,] 4
## [675,] 6
## [676,] 8
## [677,] 7
## [678,] 5
## [679,] 7
## [680,] 3
## [681,] 5
## [682,] 6
## [683,] 5
## [684,] 7
## [685,] 4
## [686,] 7
## [687,] 6
## [688,] 5
## [689,] 8
## [690,] 5
## [691,] 2
## [692,] 2
## [693,] 1
## [694,] 5
## [695,] 6
## [696,] 8
## [697,] 6
## [698,] 2
## [699,] 3
## [700,] 4
## [701,] 6
## [702,] 5
## [703,] 6
## [704,] 3
## [705,] 6
## [706,] 3
## [707,] 5
## [708,] 6
## [709,] 8
## [710,] 8
## [711,] 4
## [712,] 7
## [713,] 7
## [714,] 5
## [715,] 6
## [716,] 7
## [717,] 7
## [718,] 4
## [719,] 7
## [720,] 3
## [721,] 4
## [722,] 2
## [723,] 5
## [724,] 5
## [725,] 6
## [726,] 5
## [727,] 4
## [728,] 4
## [729,] 4
## [730,] 4
## [731,] 3
## [732,] 5
## [733,] 6
## [734,] 5
## [735,] 6
## [736,] 4
## [737,] 6
## [738,] 6
## [739,] 6
## [740,] 8
## [741,] 7
## [742,] 6
## [743,] 9
## [744,] 9
## [745,] 9
## [746,] 8
## [747,] 5
## [748,] 3
## [749,] 8
## [750,] 7
## [751,] 4
## [752,] 6
## [753,] 5
## [754,] 7
## [755,] 5
## [756,] 5
## [757,] 7
## [758,] 5
## [759,] 7
## [760,] 5
## [761,] 6
## [762,] 9
## [763,] 5
## [764,] 8
## [765,] 7
## [766,] 7
## [767,] 2
## [768,] 8
## [769,] 7
## [770,] 4
## [771,] 5
## [772,] 6
## [773,] 8
## [774,] 7
## [775,] 8
## [776,] 7
## [777,] 7
## [778,] 5
## [779,] 5
## [780,] 4
## [781,] 3
## [782,] 7
## [783,] 5
## [784,] 3
## [785,] 4
## [786,] 3
## [787,] 3
## [788,] 5
## [789,] 5
## [790,] 8
## [791,] 9
## [792,] 8
## [793,] 9
## [794,] 8
## [795,] 5
## [796,] 8
## [797,] 8
## [798,] 6
## [799,] 7
## [800,] 6
## [801,] 8
## [802,] 9
## [803,] 9
## [804,] 8
## [805,] 6
## [806,] 8
## [807,] 6
## [808,] 7
## [809,] 8
## [810,] 8
## [811,] 4
## [812,] 6
## [813,] 7
## [814,] 9
## [815,] 7
## [816,] 4
## [817,] 9
## [818,] 6
## [819,] 5
## [820,] 6
## [821,] 3
## [822,] 5
## [823,] 6
## [824,] 7
## [825,] 5
## [826,] 6
## [827,] 3
## [828,] 3
## [829,] 5
## [830,] 5
## [831,] 5
## [832,] 5
## [833,] 8
## [834,] 4
## [835,] 6
## [836,] 4
## [837,] 7
## [838,] 1
## [839,] 4
## [840,] 2
## [841,] 8
## [842,] 8
## [843,] 8
## [844,] 8
## [845,] 6
## [846,] 3
## [847,] 4
## [848,] 3
## [849,] 3
## [850,] 6
## [851,] 8
## [852,] 8
## [853,] 5
## [854,] 7
## [855,] 5
## [856,] 5
## [857,] 4
## [858,] 8
## [859,] 8
## [860,] 8
## [861,] 7
## [862,] 7
## [863,] 5
## [864,] 8
## [865,] 5
## [866,] 7
## [867,] 8
## [868,] 9
## [869,] 6
## [870,] 8
## [871,] 2
## [872,] 7
## [873,] 9
## [874,] 4
## [875,] 3
## [876,] 0
## [877,] 5
## [878,] 4
## [879,] 2
## [880,] 3
## [881,] 0
## [882,] 9
## [883,] 9
## [884,] 7
## [885,] 5
## [886,] 5
## [887,] 5
## [888,] 3
## [889,] 5
## [890,] 6
## [891,] 6
## [892,] 3
## [893,] 7
## [894,] 6
## [895,] 4
## [896,] 7
## [897,] 4
## [898,] 8
## [899,] 6
## [900,] 5
## [901,] 8
## [902,] 6
## [903,] 5
## [904,] 6
## [905,] 5
## [906,] 7
## [907,] 7
## [908,] 9
## [909,] 8
## [910,] 8
## [911,] 7
## [912,] 8
## [913,] 4
## [914,] 9
## [915,] 4
## [916,] 8
## [917,] 5
## [918,] 5
## [919,] 6
## [920,] 4
## [921,] 6
## [922,] 5
## [923,] 6
## [924,] 5
## [925,] 8
## [926,] 6
## [927,] 3
## [928,] 4
## [929,] 5
## [930,] 5
## [931,] 5
## [932,] 7
## [933,] 6
## [934,] 6
## [935,] 2
## [936,] 6
## [937,] 2
## [938,] 7
## [939,] 9
## [940,] 8
## [941,] 2
## [942,] 4
## [943,] 7
## [944,] 8
## [945,] 8
## [946,] 7
## [947,] 4
## [948,] 3
## [949,] 7
## [950,] 4
## [951,] 6
## [952,] 5
## [953,] 3
## [954,] 5
## [955,] 4
## [956,] 2
## [957,] 6
## [958,] 6
## [959,] 8
## [960,] 9
## [961,] 7
## [962,] 7
## [963,] 6
## [964,] 6
## [965,] 3
## [966,] 3
## [967,] 4
## [968,] 4
## [969,] 4
## [970,] 5
## [971,] 6
## [972,] 6
## [973,] 3
## [974,] 6
## [975,] 4
## [976,] 6
## [977,] 6
## [978,] 9
## [979,] 9
## [980,] 9
## [981,] 8
## [982,] 8
## [983,] 4
## [984,] 4
## [985,] 6
## [986,] 9
## [987,] 3
## [988,] 5
## [989,] 2
## [990,] 5
## [991,] 5
## [992,] 6
## [993,] 7
## [994,] 7
## [995,] 7
## [996,] 7
## [997,] 5
## [998,] 1
## [999,] 4
## [1000,] 3
## [1001,] 4
## [1002,] 5
## [1003,] 7
## [1004,] 9
## [1005,] 6
## [1006,] 9
## [1007,] 6
## [1008,] 8
## [1009,] 8
## [1010,] 5
## [1011,] 4
## [1012,] 6
## [1013,] 6
## [1014,] 7
## [1015,] 6
## [1016,] 4
## [1017,] 6
## [1018,] 4
## [1019,] 6
## [1020,] 3
## [1021,] 3
## [1022,] 7
## [1023,] 3
## [1024,] 7
## [1025,] 2
## [1026,] 7
## [1027,] 7
## [1028,] 8
## [1029,] 7
## [1030,] 7
## [1031,] 8
## [1032,] 7
## [1033,] 7
## [1034,] 4
## [1035,] 5
## [1036,] 5
## [1037,] 3
## [1038,] 6
## [1039,] 2
## [1040,] 7
## [1041,] 3
## [1042,] 7
## [1043,] 5
## [1044,] 3
## [1045,] 7
## [1046,] 6
## [1047,] 4
## [1048,] 6
## [1049,] 8
## [1050,] 9
## [1051,] 5
## [1052,] 6
## [1053,] 5
## [1054,] 8
## [1055,] 6
## [1056,] 7
## [1057,] 1
## [1058,] 5
## [1059,] 3
## [1060,] 7
## [1061,] 5
## [1062,] 8
## [1063,] 4
## [1064,] 9
## [1065,] 9
## [1066,] 4
## [1067,] 6
## [1068,] 7
## [1069,] 7
## [1070,] 8
## [1071,] 8
## [1072,] 4
## [1073,] 7
## [1074,] 5
## [1075,] 4
## [1076,] 7
## [1077,] 8
## [1078,] 3
## [1079,] 6
## [1080,] 9
## [1081,] 8
## [1082,] 6
## [1083,] 5
## [1084,] 8
## [1085,] 8
## [1086,] 8
## [1087,] 7
## [1088,] 8
## [1089,] 7
## [1090,] 5
## [1091,] 8
## [1092,] 6
## [1093,] 3
## [1094,] 6
## [1095,] 7
## [1096,] 6
## [1097,] 7
## [1098,] 8
## [1099,] 3
## [1100,] 5
## [1101,] 3
## [1102,] 8
## [1103,] 1
## [1104,] 5
## [1105,] 8
## [1106,] 7
## [1107,] 8
## [1108,] 3
## [1109,] 3
## [1110,] 7
## [1111,] 8
## [1112,] 9
## [1113,] 7
## [1114,] 7
## [1115,] 9
## [1116,] 8
## [1117,] 8
## [1118,] 6
## [1119,] 8
## [1120,] 7
## [1121,] 5
## [1122,] 7
## [1123,] 3
## [1124,] 8
## [1125,] 8
## [1126,] 5
## [1127,] 5
## [1128,] 6
## [1129,] 4
## [1130,] 3
## [1131,] 5
## [1132,] 6
## [1133,] 8
## [1134,] 5
## [1135,] 7
## [1136,] 5
## [1137,] 5
## [1138,] 6
## [1139,] 2
## [1140,] 2
## [1141,] 7
## [1142,] 4
## [1143,] 8
## [1144,] 8
## [1145,] 3
## [1146,] 9
## [1147,] 9
## [1148,] 7
## [1149,] 7
## [1150,] 6
## [1151,] 6
## [1152,] 5
## [1153,] 6
## [1154,] 9
## [1155,] 9
## [1156,] 9
## [1157,] 7
## [1158,] 2
## [1159,] 8
## [1160,] 6
## [1161,] 4
## [1162,] 7
## [1163,] 2
## [1164,] 6
## [1165,] 6
## [1166,] 3
## [1167,] 5
## [1168,] 3
## [1169,] 6
## [1170,] 8
## [1171,] 6
## [1172,] 7
## [1173,] 3
## [1174,] 2
## [1175,] 4
## [1176,] 9
## [1177,] 8
## [1178,] 9
## [1179,] 7
## [1180,] 8
## [1181,] 4
## [1182,] 6
## [1183,] 6
## [1184,] 9
## [1185,] 8
## [1186,] 8
## [1187,] 4
## [1188,] 4
## [1189,] 6
## [1190,] 6
## [1191,] 3
## [1192,] 7
## [1193,] 3
## [1194,] 5
## [1195,] 7
## [1196,] 5
## [1197,] 7
## [1198,] 8
## [1199,] 9
## [1200,] 5
## [1201,] 5
## [1202,] 9
## [1203,] 2
## [1204,] 2
## [1205,] 4
## [1206,] 3
## [1207,] 4
## [1208,] 5
## [1209,] 7
## [1210,] 8
## [1211,] 6
## [1212,] 7
## [1213,] 8
## [1214,] 8
## [1215,] 8
## [1216,] 7
## [1217,] 5
## [1218,] 7
## [1219,] 8
## [1220,] 8
## [1221,] 6
## [1222,] 8
## [1223,] 5
## [1224,] 7
## [1225,] 6
## [1226,] 6
## [1227,] 6
## [1228,] 5
## [1229,] 4
## [1230,] 6
## [1231,] 5
## [1232,] 7
## [1233,] 5
## [1234,] 1
## [1235,] 3
## [1236,] 3
## [1237,] 3
## [1238,] 2
## [1239,] 8
## [1240,] 6
## [1241,] 9
## [1242,] 9
## [1243,] 9
## [1244,] 8
## [1245,] 6
## [1246,] 9
## [1247,] 8
## [1248,] 7
## [1249,] 5
## [1250,] 5
## [1251,] 7
## [1252,] 5
## [1253,] 9
## [1254,] 7
## [1255,] 9
## [1256,] 7
## [1257,] 6
## [1258,] 5
## [1259,] 7
## [1260,] 7
## [1261,] 9
## [1262,] 6
## [1263,] 2
## [1264,] 3
## [1265,] 2
## [1266,] 3
## [1267,] 3
## [1268,] 7
## [1269,] 6
## [1270,] 6
## [1271,] 5
## [1272,] 6
## [1273,] 6
## [1274,] 4
## [1275,] 7
## [1276,] 7
## [1277,] 5
## [1278,] 4
## [1279,] 6
## [1280,] 7
## [1281,] 4
## [1282,] 8
## [1283,] 3
## [1284,] 5
## [1285,] 2
## [1286,] 4
## [1287,] 6
## [1288,] 3
## [1289,] 5
## [1290,] 8
## [1291,] 6
## [1292,] 9
## [1293,] 5
## [1294,] 7
## [1295,] 6
## [1296,] 8
## [1297,] 4
## [1298,] 8
## [1299,] 8
## [1300,] 6
## [1301,] 6
## [1302,] 5
## [1303,] 7
## [1304,] 5
## [1305,] 7
## [1306,] 9
## [1307,] 7
## [1308,] 8
## [1309,] 7
## [1310,] 8
## [1311,] 6
## [1312,] 9
## [1313,] 8
## [1314,] 5
## [1315,] 4
## [1316,] 4
## [1317,] 7
## [1318,] 4
## [1319,] 8
## [1320,] 7
## [1321,] 6
## [1322,] 7
## [1323,] 5
## [1324,] 6
## [1325,] 6
## [1326,] 7
## [1327,] 4
## [1328,] 8
## [1329,] 8
## [1330,] 6
## [1331,] 8
## [1332,] 6
## [1333,] 4
## [1334,] 5
## [1335,] 4
## [1336,] 5
## [1337,] 6
## [1338,] 3
## [1339,] 7
## [1340,] 8
## [1341,] 8
## [1342,] 9
## [1343,] 8
## [1344,] 7
## [1345,] 6
## [1346,] 6
## [1347,] 7
## [1348,] 4
## [1349,] 6
## [1350,] 4
## [1351,] 8
## [1352,] 4
## [1353,] 7
## [1354,] 6
## [1355,] 7
## [1356,] 7
## [1357,] 6
## [1358,] 3
## [1359,] 5
## [1360,] 2
## [1361,] 5
## [1362,] 3
## [1363,] 2
## [1364,] 6
## [1365,] 2
## [1366,] 6
## [1367,] 6
## [1368,] 3
## [1369,] 6
## [1370,] 5
## [1371,] 5
## [1372,] 1
## [1373,] 7
## [1374,] 4
## [1375,] 8
## [1376,] 6
## [1377,] 7
## [1378,] 2
## [1379,] 3
## [1380,] 3
## [1381,] 5
## [1382,] 6
## [1383,] 4
## [1384,] 5
## [1385,] 8
## [1386,] 7
## [1387,] 8
## [1388,] 3
## [1389,] 7
## [1390,] 6
## [1391,] 9
## [1392,] 4
## [1393,] 7
## [1394,] 4
## [1395,] 4
## [1396,] 8
## [1397,] 7
## [1398,] 6
## [1399,] 6
## [1400,] 5
## [1401,] 7
## [1402,] 3
## [1403,] 7
## [1404,] 6
## [1405,] 7
## [1406,] 5
## [1407,] 6
## [1408,] 6
## [1409,] 3
## [1410,] 9
## [1411,] 9
## [1412,] 7
## [1413,] 4
## [1414,] 6
## [1415,] 5
## [1416,] 6
## [1417,] 7
## [1418,] 8
## [1419,] 6
## [1420,] 8
## [1421,] 8
## [1422,] 8
## [1423,] 4
## [1424,] 6
## [1425,] 5
## [1426,] 6
## [1427,] 4
## [1428,] 5
## [1429,] 5
## [1430,] 6
## [1431,] 4
## [1432,] 3
## [1433,] 4
## [1434,] 4
## [1435,] 5
## [1436,] 3
## [1437,] 2
## [1438,] 7
## [1439,] 3
## [1440,] 6
## [1441,] 6
## [1442,] 4
## [1443,] 6
## [1444,] 7
## [1445,] 4
## [1446,] 6
## [1447,] 6
## [1448,] 8
## [1449,] 7
## [1450,] 2
## [1451,] 7
## [1452,] 7
## [1453,] 7
## [1454,] 3
## [1455,] 8
## [1456,] 7
## [1457,] 8
## [1458,] 5
## [1459,] 8
## [1460,] 3
## [1461,] 4
## [1462,] 6
## [1463,] 2
## [1464,] 5
## [1465,] 3
## [1466,] 4
## [1467,] 7
## [1468,] 6
## [1469,] 4
## [1470,] 6
## [1471,] 4
## [1472,] 9
## [1473,] 5
## [1474,] 7
## [1475,] 3
## [1476,] 4
## [1477,] 7
## [1478,] 5
## [1479,] 9
## [1480,] 5
## [1481,] 7
## [1482,] 8
## [1483,] 5
## [1484,] 8
## [1485,] 8
## [1486,] 9
## [1487,] 3
## [1488,] 3
## [1489,] 7
## [1490,] 6
## [1491,] 5
## [1492,] 5
## [1493,] 5
## [1494,] 6
## [1495,] 6
## [1496,] 5
## [1497,] 7
## [1498,] 6
## [1499,] 5
## [1500,] 6
## [1501,] 6
## [1502,] 7
## [1503,] 7
## [1504,] 5
## [1505,] 7
## [1506,] 8
## [1507,] 7
## [1508,] 7
## [1509,] 6
## [1510,] 4
## [1511,] 4
## [1512,] 5
## [1513,] 2
## [1514,] 6
## [1515,] 8
## [1516,] 7
## [1517,] 2
## [1518,] 5
## [1519,] 5
## [1520,] 4
## [1521,] 3
## [1522,] 1
## [1523,] 7
## [1524,] 6
## [1525,] 7
## [1526,] 6
## [1527,] 5
## [1528,] 5
## [1529,] 9
## [1530,] 8
## [1531,] 4
## [1532,] 5
## [1533,] 5
## [1534,] 6
## [1535,] 7
## [1536,] 4
## [1537,] 8
## [1538,] 5
## [1539,] 5
## [1540,] 5
## [1541,] 5
## [1542,] 4
## [1543,] 4
## [1544,] 3
## [1545,] 2
## [1546,] 5
## [1547,] 3
## [1548,] 6
## [1549,] 6
## [1550,] 7
## [1551,] 8
## [1552,] 3
## [1553,] 4
## [1554,] 5
## [1555,] 6
## [1556,] 5
## [1557,] 6
## [1558,] 5
## [1559,] 5
## [1560,] 6
## [1561,] 8
## [1562,] 5
## [1563,] 6
## [1564,] 4
## [1565,] 5
## [1566,] 6
## [1567,] 6
## [1568,] 7
## [1569,] 8
## [1570,] 5
## [1571,] 2
## [1572,] 5
## [1573,] 9
## [1574,] 9
## [1575,] 8
## [1576,] 8
## [1577,] 7
## [1578,] 5
## [1579,] 2
## [1580,] 2
## [1581,] 5
## [1582,] 4
## [1583,] 5
## [1584,] 4
## [1585,] 6
## [1586,] 7
## [1587,] 5
## [1588,] 7
## [1589,] 5
## [1590,] 6
## [1591,] 5
## [1592,] 2
## [1593,] 7
## [1594,] 5
## [1595,] 5
## [1596,] 5
## [1597,] 7
## [1598,] 7
## [1599,] 6
## [1600,] 6
## [1601,] 7
## [1602,] 6
## [1603,] 6
## [1604,] 8
## [1605,] 7
## [1606,] 7
## [1607,] 6
## [1608,] 6
## [1609,] 8
## [1610,] 6
## [1611,] 6
## [1612,] 2
## [1613,] 5
## [1614,] 5
## [1615,] 5
## [1616,] 4
## [1617,] 4
## [1618,] 7
## [1619,] 4
## [1620,] 6
## [1621,] 7
## [1622,] 7
## [1623,] 6
## [1624,] 9
## [1625,] 6
## [1626,] 2
## [1627,] 2
## [1628,] 5
## [1629,] 4
## [1630,] 8
## [1631,] 6
## [1632,] 4
## [1633,] 5
## [1634,] 8
## [1635,] 3
## [1636,] 5
## [1637,] 7
## [1638,] 3
## [1639,] 4
## [1640,] 4
## [1641,] 4
## [1642,] 7
## [1643,] 5
## [1644,] 4
## [1645,] 6
## [1646,] 7
## [1647,] 2
## [1648,] 5
## [1649,] 7
## [1650,] 8
## [1651,] 7
## [1652,] 4
## [1653,] 5
## [1654,] 8
## [1655,] 8
## [1656,] 5
## [1657,] 6
## [1658,] 6
## [1659,] 6
## [1660,] 4
## [1661,] 5
## [1662,] 1
## [1663,] 4
## [1664,] 3
## [1665,] 2
## [1666,] 6
## [1667,] 3
## [1668,] 2
## [1669,] 8
## [1670,] 8
## [1671,] 4
## [1672,] 8
## [1673,] 5
## [1674,] 6
## [1675,] 5
## [1676,] 5
## [1677,] 7
## [1678,] 6
## [1679,] 2
## [1680,] 1
## [1681,] 6
## [1682,] 5
## [1683,] 6
## [1684,] 6
## [1685,] 7
## [1686,] 7
## [1687,] 4
## [1688,] 9
## [1689,] 4
## [1690,] 5
## [1691,] 5
## [1692,] 5
## [1693,] 2
## [1694,] 5
## [1695,] 8
## [1696,] 6
## [1697,] 8
## [1698,] 4
## [1699,] 8
## [1700,] 8
## [1701,] 5
## [1702,] 5
## [1703,] 9
## [1704,] 7
## [1705,] 7
## [1706,] 2
## [1707,] 8
## [1708,] 5
## [1709,] 7
## [1710,] 6
## [1711,] 8
## [1712,] 5
## [1713,] 5
## [1714,] 4
## [1715,] 8
## [1716,] 5
## [1717,] 3
## [1718,] 4
## [1719,] 6
## [1720,] 8
## [1721,] 7
## [1722,] 8
## [1723,] 7
## [1724,] 5
## [1725,] 5
## [1726,] 5
## [1727,] 7
## [1728,] 8
## [1729,] 7
## [1730,] 5
## [1731,] 5
## [1732,] 1
## [1733,] 5
## [1734,] 5
## [1735,] 5
## [1736,] 8
## [1737,] 6
## [1738,] 5
## [1739,] 5
## [1740,] 4
## [1741,] 9
## [1742,] 6
## [1743,] 7
## [1744,] 5
## [1745,] 7
## [1746,] 7
## [1747,] 8
## [1748,] 7
## [1749,] 9
## [1750,] 8
## [1751,] 5
## [1752,] 5
## [1753,] 7
## [1754,] 5
## [1755,] 5
## [1756,] 7
## [1757,] 9
## [1758,] 2
## [1759,] 2
## [1760,] 3
## [1761,] 2
## [1762,] 7
## [1763,] 6
## [1764,] 6
## [1765,] 7
## [1766,] 6
## [1767,] 7
## [1768,] 2
## [1769,] 4
## [1770,] 6
## [1771,] 7
## [1772,] 4
## [1773,] 4
## [1774,] 2
## [1775,] 5
## [1776,] 6
## [1777,] 6
## [1778,] 7
## [1779,] 7
## [1780,] 5
## [1781,] 8
## [1782,] 8
## [1783,] 8
## [1784,] 9
## [1785,] 9
## [1786,] 7
## [1787,] 6
## [1788,] 5
## [1789,] 5
## [1790,] 6
## [1791,] 5
## [1792,] 6
## [1793,] 5
## [1794,] 5
## [1795,] 8
## [1796,] 6
## [1797,] 9
## [1798,] 6
## [1799,] 4
## [1800,] 5
## [1801,] 3
## [1802,] 6
## [1803,] 8
## [1804,] 7
## [1805,] 8
## [1806,] 7
## [1807,] 5
## [1808,] 9
## [1809,] 6
## [1810,] 7
## [1811,] 8
## [1812,] 7
## [1813,] 3
## [1814,] 5
## [1815,] 6
## [1816,] 5
## [1817,] 2
## [1818,] 8
## [1819,] 2
## [1820,] 3
## [1821,] 6
## [1822,] 3
## [1823,] 2
## [1824,] 7
## [1825,] 4
## [1826,] 8
## [1827,] 6
## [1828,] 6
## [1829,] 7
## [1830,] 7
## [1831,] 8
## [1832,] 7
## [1833,] 7
## [1834,] 2
## [1835,] 8
## [1836,] 3
## [1837,] 2
## [1838,] 5
## [1839,] 8
## [1840,] 5
## [1841,] 2
## [1842,] 5
## [1843,] 7
## [1844,] 2
## [1845,] 7
## [1846,] 4
## [1847,] 8
## [1848,] 7
## [1849,] 8
## [1850,] 5
## [1851,] 7
## [1852,] 7
## [1853,] 6
## [1854,] 2
## [1855,] 7
## [1856,] 7
## [1857,] 4
## [1858,] 3
## [1859,] 6
## [1860,] 7
## [1861,] 7
## [1862,] 5
## [1863,] 2
## [1864,] 4
## [1865,] 7
## [1866,] 4
## [1867,] 4
## [1868,] 8
## [1869,] 8
## [1870,] 8
## [1871,] 7
## [1872,] 7
## [1873,] 5
## [1874,] 9
## [1875,] 8
## [1876,] 7
## [1877,] 8
## [1878,] 3
## [1879,] 5
## [1880,] 4
## [1881,] 8
## [1882,] 9
## [1883,] 9
## [1884,] 7
## [1885,] 3
## [1886,] 6
## [1887,] 4
## [1888,] 7
## [1889,] 8
## [1890,] 5
## [1891,] 5
## [1892,] 7
## [1893,] 9
## [1894,] 5
## [1895,] 7
## [1896,] 6
## [1897,] 2
## [1898,] 7
## [1899,] 6
## [1900,] 6
## [1901,] 5
## [1902,] 5
## [1903,] 7
## [1904,] 8
## [1905,] 5
## [1906,] 7
## [1907,] 4
## [1908,] 6
## [1909,] 7
## [1910,] 9
## [1911,] 8
## [1912,] 7
## [1913,] 8
## [1914,] 5
## [1915,] 6
## [1916,] 7
## [1917,] 7
## [1918,] 5
## [1919,] 5
## [1920,] 8
## [1921,] 4
## [1922,] 6
## [1923,] 5
## [1924,] 4
## [1925,] 9
## [1926,] 7
## [1927,] 4
## [1928,] 4
## [1929,] 9
## [1930,] 5
## [1931,] 5
## [1932,] 6
## [1933,] 8
## [1934,] 5
## [1935,] 6
## [1936,] 7
## [1937,] 2
## [1938,] 2
## [1939,] 3
## [1940,] 4
## [1941,] 3
## [1942,] 3
## [1943,] 4
## [1944,] 3
## [1945,] 8
## [1946,] 5
## [1947,] 4
## [1948,] 3
## [1949,] 5
## [1950,] 5
## [1951,] 8
## [1952,] 8
## [1953,] 5
## [1954,] 4
## [1955,] 3
## [1956,] 7
## [1957,] 8
## [1958,] 7
## [1959,] 5
## [1960,] 9
## [1961,] 8
## [1962,] 1
## [1963,] 8
## [1964,] 4
## [1965,] 7
## [1966,] 6
## [1967,] 6
## [1968,] 6
## [1969,] 6
## [1970,] 7
## [1971,] 5
## [1972,] 6
## [1973,] 3
## [1974,] 5
## [1975,] 7
## [1976,] 7
## [1977,] 5
## [1978,] 5
## [1979,] 5
## [1980,] 5
## [1981,] 6
## [1982,] 3
## [1983,] 5
## [1984,] 5
## [1985,] 7
## [1986,] 8
## [1987,] 7
## [1988,] 8
## [1989,] 7
## [1990,] 9
## [1991,] 5
## [1992,] 5
## [1993,] 6
## [1994,] 5
## [1995,] 1
## [1996,] 3
## [1997,] 2
## [1998,] 2
## [1999,] 4
## [2000,] 9
## [2001,] 8
## [2002,] 4
## [2003,] 6
## [2004,] 6
## [2005,] 5
## [2006,] 6
## [2007,] 8
## [2008,] 9
## [2009,] 7
## [2010,] 6
## [2011,] 3
## [2012,] 6
## [2013,] 7
## [2014,] 8
## [2015,] 6
## [2016,] 9
## [2017,] 6
## [2018,] 5
## [2019,] 6
## [2020,] 5
## [2021,] 8
## [2022,] 2
## [2023,] 8
## [2024,] 9
## [2025,] 5
## [2026,] 7
## [2027,] 6
## [2028,] 2
## [2029,] 7
## [2030,] 6
## [2031,] 5
## [2032,] 9
## [2033,] 6
## [2034,] 5
## [2035,] 5
## [2036,] 4
## [2037,] 5
## [2038,] 5
## [2039,] 6
## [2040,] 5
## [2041,] 5
## [2042,] 3
## [2043,] 6
## [2044,] 7
## [2045,] 4
## [2046,] 7
## [2047,] 7
## [2048,] 5
## [2049,] 2
## [2050,] 2
## [2051,] 5
## [2052,] 3
## [2053,] 5
## [2054,] 2
## [2055,] 3
## [2056,] 8
## [2057,] 6
## [2058,] 5
## [2059,] 5
## [2060,] 4
## [2061,] 8
## [2062,] 3
## [2063,] 2
## [2064,] 8
## [2065,] 8
## [2066,] 8
## [2067,] 6
## [2068,] 6
## [2069,] 5
## [2070,] 9
## [2071,] 8
## [2072,] 8
## [2073,] 5
## [2074,] 4
## [2075,] 4
## [2076,] 2
## [2077,] 4
## [2078,] 3
## [2079,] 3
## [2080,] 4
## [2081,] 5
## [2082,] 3
## [2083,] 5
## [2084,] 6
## [2085,] 7
## [2086,] 5
## [2087,] 6
## [2088,] 6
## [2089,] 7
## [2090,] 5
## [2091,] 6
## [2092,] 8
## [2093,] 6
## [2094,] 4
## [2095,] 4
## [2096,] 5
## [2097,] 5
## [2098,] 9
## [2099,] 9
## [2100,] 7
## [2101,] 4
## [2102,] 7
## [2103,] 5
## [2104,] 3
## [2105,] 5
## [2106,] 7
## [2107,] 7
## [2108,] 8
## [2109,] 8
## [2110,] 7
## [2111,] 7
## [2112,] 7
## [2113,] 7
## [2114,] 8
## [2115,] 7
## [2116,] 6
## [2117,] 8
## [2118,] 8
## [2119,] 6
## [2120,] 4
## [2121,] 5
## [2122,] 4
## [2123,] 2
## [2124,] 1
## [2125,] 5
## [2126,] 4
## [2127,] 6
## [2128,] 7
## [2129,] 6
## [2130,] 3
## [2131,] 4
## [2132,] 3
## [2133,] 5
## [2134,] 6
## [2135,] 5
## [2136,] 5
## [2137,] 7
## [2138,] 7
## [2139,] 7
## [2140,] 6
## [2141,] 7
## [2142,] 6
## [2143,] 7
## [2144,] 4
## [2145,] 7
## [2146,] 7
## [2147,] 3
## [2148,] 5
## [2149,] 6
## [2150,] 5
## [2151,] 8
## [2152,] 6
## [2153,] 4
## [2154,] 7
## [2155,] 5
## [2156,] 8
## [2157,] 6
## [2158,] 5
## [2159,] 6
## [2160,] 8
## [2161,] 4
## [2162,] 6
## [2163,] 6
## [2164,] 9
## [2165,] 5
## [2166,] 4
## [2167,] 5
## [2168,] 4
## [2169,] 7
## [2170,] 7
## [2171,] 6
## [2172,] 2
## [2173,] 8
## [2174,] 6
## [2175,] 7
## [2176,] 4
## [2177,] 5
## [2178,] 6
## [2179,] 7
## [2180,] 4
## [2181,] 7
## [2182,] 5
## [2183,] 3
## [2184,] 5
## [2185,] 3
## [2186,] 3
## [2187,] 5
## [2188,] 8
## [2189,] 9
## [2190,] 9
## [2191,] 7
## [2192,] 6
## [2193,] 7
## [2194,] 5
## [2195,] 8
## [2196,] 6
## [2197,] 3
## [2198,] 9
## [2199,] 9
## [2200,] 8
## [2201,] 7
## [2202,] 5
## [2203,] 6
## [2204,] 6
## [2205,] 4
## [2206,] 5
## [2207,] 5
## [2208,] 6
## [2209,] 5
## [2210,] 6
## [2211,] 4
## [2212,] 8
## [2213,] 4
## [2214,] 3
## [2215,] 4
## [2216,] 2
## [2217,] 5
## [2218,] 7
## [2219,] 4
## [2220,] 6
## [2221,] 6
## [2222,] 4
## [2223,] 5
## [2224,] 3
## [2225,] 4
## [2226,] 5
## [2227,] 5
## [2228,] 6
## [2229,] 4
## [2230,] 4
## [2231,] 5
## [2232,] 7
## [2233,] 7
## [2234,] 6
## [2235,] 5
## [2236,] 5
## [2237,] 4
## [2238,] 6
## [2239,] 5
## [2240,] 8
## [2241,] 9
## [2242,] 7
## [2243,] 6
## [2244,] 7
## [2245,] 4
## [2246,] 6
## [2247,] 6
## [2248,] 8
## [2249,] 5
## [2250,] 9
## [2251,] 7
## [2252,] 5
## [2253,] 8
## [2254,] 4
## [2255,] 4
## [2256,] 3
## [2257,] 5
## [2258,] 8
## [2259,] 5
## [2260,] 8
## [2261,] 5
## [2262,] 6
## [2263,] 7
## [2264,] 6
## [2265,] 7
## [2266,] 8
## [2267,] 5
## [2268,] 4
## [2269,] 5
## [2270,] 4
## [2271,] 8
## [2272,] 6
## [2273,] 7
## [2274,] 4
## [2275,] 7
## [2276,] 9
## [2277,] 6
## [2278,] 6
## [2279,] 4
## [2280,] 4
## [2281,] 5
## [2282,] 7
## [2283,] 9
## [2284,] 4
## [2285,] 5
## [2286,] 6
## [2287,] 9
## [2288,] 7
## [2289,] 8
## [2290,] 5
## [2291,] 7
## [2292,] 8
## [2293,] 6
## [2294,] 8
## [2295,] 8
## [2296,] 7
## [2297,] 5
## [2298,] 7
## [2299,] 6
## [2300,] 8
## [2301,] 5
## [2302,] 6
## [2303,] 7
## [2304,] 5
## [2305,] 5
## [2306,] 5
## [2307,] 4
## [2308,] 3
## [2309,] 6
## [2310,] 6
## [2311,] 5
## [2312,] 7
## [2313,] 5
## [2314,] 5
## [2315,] 2
## [2316,] 4
## [2317,] 3
## [2318,] 5
## [2319,] 4
## [2320,] 6
## [2321,] 5
## [2322,] 5
## [2323,] 5
## [2324,] 6
## [2325,] 7
## [2326,] 7
## [2327,] 6
## [2328,] 6
## [2329,] 9
## [2330,] 5
## [2331,] 7
## [2332,] 4
## [2333,] 6
## [2334,] 6
## [2335,] 8
## [2336,] 6
## [2337,] 6
## [2338,] 6
## [2339,] 6
## [2340,] 4
## [2341,] 8
## [2342,] 9
## [2343,] 9
## [2344,] 9
## [2345,] 6
## [2346,] 3
## [2347,] 8
## [2348,] 8
## [2349,] 5
## [2350,] 6
## [2351,] 4
## [2352,] 8
## [2353,] 5
## [2354,] 4
## [2355,] 4
## [2356,] 7
## [2357,] 7
## [2358,] 5
## [2359,] 5
## [2360,] 3
## [2361,] 6
## [2362,] 8
## [2363,] 8
## [2364,] 7
## [2365,] 4
## [2366,] 1
## [2367,] 3
## [2368,] 5
## [2369,] 5
## [2370,] 6
## [2371,] 8
## [2372,] 9
## [2373,] 8
## [2374,] 1
## [2375,] 6
## [2376,] 6
## [2377,] 7
## [2378,] 6
## [2379,] 5
## [2380,] 4
## [2381,] 6
## [2382,] 8
## [2383,] 9
## [2384,] 8
## [2385,] 5
## [2386,] 6
## [2387,] 3
## [2388,] 6
## [2389,] 1
## [2390,] 4
## [2391,] 8
## [2392,] 7
## [2393,] 4
## [2394,] 3
## [2395,] 9
## [2396,] 7
## [2397,] 5
## [2398,] 7
## [2399,] 4
## [2400,] 5
## [2401,] 8
## [2402,] 8
## [2403,] 8
## [2404,] 7
## [2405,] 8
## [2406,] 5
## [2407,] 8
## [2408,] 6
## [2409,] 5
## [2410,] 3
## [2411,] 5
## [2412,] 7
## [2413,] 7
## [2414,] 3
## [2415,] 6
## [2416,] 6
## [2417,] 7
## [2418,] 6
## [2419,] 8
## [2420,] 7
## [2421,] 4
## [2422,] 4
## [2423,] 6
## [2424,] 4
## [2425,] 6
## [2426,] 3
## [2427,] 5
## [2428,] 6
## [2429,] 7
## [2430,] 5
## [2431,] 7
## [2432,] 3
## [2433,] 7
## [2434,] 4
## [2435,] 6
## [2436,] 8
## [2437,] 5
## [2438,] 5
## [2439,] 5
## [2440,] 5
## [2441,] 3
## [2442,] 4
## [2443,] 4
## [2444,] 3
## [2445,] 4
## [2446,] 3
## [2447,] 4
## [2448,] 6
## [2449,] 5
## [2450,] 4
## [2451,] 4
## [2452,] 7
## [2453,] 3
## [2454,] 4
## [2455,] 9
## [2456,] 6
## [2457,] 9
## [2458,] 4
## [2459,] 7
## [2460,] 5
## [2461,] 5
## [2462,] 7
## [2463,] 7
## [2464,] 3
## [2465,] 3
## [2466,] 4
## [2467,] 4
## [2468,] 8
## [2469,] 8
## [2470,] 3
## [2471,] 6
## [2472,] 5
## [2473,] 6
## [2474,] 5
## [2475,] 7
## [2476,] 8
## [2477,] 8
## [2478,] 7
## [2479,] 7
## [2480,] 5
## [2481,] 6
## [2482,] 8
## [2483,] 7
## [2484,] 6
## [2485,] 4
## [2486,] 5
## [2487,] 5
## [2488,] 5
## [2489,] 6
## [2490,] 7
## [2491,] 3
## [2492,] 7
## [2493,] 4
## [2494,] 6
## [2495,] 7
## [2496,] 7
## [2497,] 3
## [2498,] 6
## [2499,] 7
## [2500,] 6
## [2501,] 7
## [2502,] 6
## [2503,] 6
## [2504,] 7
## [2505,] 8
## [2506,] 6
## [2507,] 7
## [2508,] 9
## [2509,] 5
## [2510,] 7
## [2511,] 8
## [2512,] 3
## [2513,] 6
## [2514,] 9
## [2515,] 8
## [2516,] 7
## [2517,] 6
## [2518,] 8
## [2519,] 8
## [2520,] 7
## [2521,] 6
## [2522,] 6
## [2523,] 6
## [2524,] 8
## [2525,] 6
## [2526,] 7
## [2527,] 4
## [2528,] 7
## [2529,] 6
## [2530,] 8
## [2531,] 7
## [2532,] 7
## [2533,] 9
## [2534,] 7
## [2535,] 3
## [2536,] 6
## [2537,] 6
## [2538,] 4
## [2539,] 8
## [2540,] 8
## [2541,] 6
## [2542,] 5
## [2543,] 2
## [2544,] 6
## [2545,] 3
## [2546,] 3
## [2547,] 5
## [2548,] 8
## [2549,] 5
## [2550,] 4
## [2551,] 5
## [2552,] 6
## [2553,] 6
## [2554,] 5
## [2555,] 4
## [2556,] 8
## [2557,] 6
## [2558,] 8
## [2559,] 8
## [2560,] 3
## [2561,] 8
## [2562,] 8
## [2563,] 4
## [2564,] 8
## [2565,] 7
## [2566,] 7
## [2567,] 8
## [2568,] 8
## [2569,] 6
## [2570,] 6
## [2571,] 8
## [2572,] 5
## [2573,] 7
## [2574,] 6
## [2575,] 3
## [2576,] 7
## [2577,] 5
## [2578,] 7
## [2579,] 9
## [2580,] 3
## [2581,] 6
## [2582,] 4
## [2583,] 9
## [2584,] 4
## [2585,] 6
## [2586,] 7
## [2587,] 8
## [2588,] 7
## [2589,] 8
## [2590,] 6
## [2591,] 8
## [2592,] 4
## [2593,] 6
## [2594,] 3
## [2595,] 5
## [2596,] 6
## [2597,] 5
## [2598,] 6
## [2599,] 5
## [2600,] 2
## [2601,] 4
## [2602,] 5
## [2603,] 4
## [2604,] 6
## [2605,] 5
## [2606,] 5
## [2607,] 5
## [2608,] 5
## [2609,] 9
## [2610,] 4
## [2611,] 7
## [2612,] 5
## [2613,] 1
## [2614,] 7
## [2615,] 6
## [2616,] 5
## [2617,] 3
## [2618,] 5
## [2619,] 6
## [2620,] 7
## [2621,] 3
## [2622,] 5
## [2623,] 5
## [2624,] 4
## [2625,] 5
## [2626,] 6
## [2627,] 5
## [2628,] 6
## [2629,] 8
## [2630,] 8
## [2631,] 5
## [2632,] 4
## [2633,] 4
## [2634,] 4
## [2635,] 4
## [2636,] 6
## [2637,] 7
## [2638,] 4
## [2639,] 7
## [2640,] 6
## [2641,] 7
## [2642,] 8
## [2643,] 8
## [2644,] 8
## [2645,] 7
## [2646,] 7
## [2647,] 4
## [2648,] 4
## [2649,] 9
## [2650,] 7
## [2651,] 3
## [2652,] 2
## [2653,] 2
## [2654,] 3
## [2655,] 4
## [2656,] 6
## [2657,] 5
## [2658,] 2
## [2659,] 2
## [2660,] 7
## [2661,] 5
## [2662,] 7
## [2663,] 5
## [2664,] 6
## [2665,] 6
## [2666,] 6
## [2667,] 8
## [2668,] 7
## [2669,] 8
## [2670,] 8
## [2671,] 4
## [2672,] 8
## [2673,] 4
## [2674,] 5
## [2675,] 7
## [2676,] 7
## [2677,] 6
## [2678,] 5
## [2679,] 7
## [2680,] 6
## [2681,] 6
## [2682,] 8
## [2683,] 4
## [2684,] 7
## [2685,] 7
## [2686,] 6
## [2687,] 6
## [2688,] 5
## [2689,] 9
## [2690,] 7
## [2691,] 5
## [2692,] 7
## [2693,] 7
## [2694,] 8
## [2695,] 5
## [2696,] 6
## [2697,] 3
## [2698,] 1
## [2699,] 5
## [2700,] 4
## [2701,] 5
## [2702,] 2
## [2703,] 4
## [2704,] 5
## [2705,] 2
## [2706,] 1
## [2707,] 3
## [2708,] 6
## [2709,] 4
## [2710,] 8
## [2711,] 5
## [2712,] 4
## [2713,] 5
## [2714,] 5
## [2715,] 8
## [2716,] 6
## [2717,] 6
## [2718,] 5
## [2719,] 4
## [2720,] 5
## [2721,] 7
## [2722,] 3
## [2723,] 3
## [2724,] 4
## [2725,] 8
## [2726,] 8
## [2727,] 6
## [2728,] 7
## [2729,] 5
## [2730,] 5
## [2731,] 6
## [2732,] 6
## [2733,] 7
## [2734,] 7
## [2735,] 5
## [2736,] 7
## [2737,] 5
## [2738,] 5
## [2739,] 7
## [2740,] 6
## [2741,] 4
## [2742,] 4
## [2743,] 3
## [2744,] 4
## [2745,] 2
## [2746,] 8
## [2747,] 6
## [2748,] 6
## [2749,] 8
## [2750,] 6
## [2751,] 8
## [2752,] 5
## [2753,] 3
## [2754,] 5
## [2755,] 5
## [2756,] 7
## [2757,] 9
## [2758,] 5
## [2759,] 7
## [2760,] 5
## [2761,] 8
## [2762,] 9
## [2763,] 8
## [2764,] 7
## [2765,] 8
## [2766,] 6
## [2767,] 5
## [2768,] 6
## [2769,] 4
## [2770,] 4
## [2771,] 6
## [2772,] 5
## [2773,] 7
## [2774,] 5
## [2775,] 3
## [2776,] 9
## [2777,] 8
## [2778,] 5
## [2779,] 8
## [2780,] 6
## [2781,] 8
## [2782,] 5
## [2783,] 6
## [2784,] 6
## [2785,] 3
## [2786,] 6
## [2787,] 5
## [2788,] 6
## [2789,] 4
## [2790,] 8
## [2791,] 4
## [2792,] 3
## [2793,] 5
## [2794,] 1
## [2795,] 4
## [2796,] 3
## [2797,] 7
## [2798,] 6
## [2799,] 6
## [2800,] 7
## [2801,] 5
## [2802,] 6
## [2803,] 5
## [2804,] 3
## [2805,] 5
## [2806,] 7
## [2807,] 8
## [2808,] 7
## [2809,] 6
## [2810,] 5
## [2811,] 2
## [2812,] 5
## [2813,] 6
## [2814,] 6
## [2815,] 9
## [2816,] 7
## [2817,] 8
## [2818,] 7
## [2819,] 5
## [2820,] 7
## [2821,] 8
## [2822,] 8
## [2823,] 7
## [2824,] 5
## [2825,] 6
## [2826,] 6
## [2827,] 1
## [2828,] 4
## [2829,] 3
## [2830,] 4
## [2831,] 7
## [2832,] 7
## [2833,] 7
## [2834,] 6
## [2835,] 7
## [2836,] 8
## [2837,] 7
## [2838,] 6
## [2839,] 6
## [2840,] 5
## [2841,] 5
## [2842,] 5
## [2843,] 6
## [2844,] 6
## [2845,] 6
## [2846,] 5
## [2847,] 4
## [2848,] 4
## [2849,] 5
## [2850,] 6
## [2851,] 5
## [2852,] 6
## [2853,] 5
## [2854,] 7
## [2855,] 4
## [2856,] 5
## [2857,] 4
## [2858,] 2
## [2859,] 3
## [2860,] 6
## [2861,] 5
## [2862,] 3
## [2863,] 5
## [2864,] 3
## [2865,] 6
## [2866,] 3
## [2867,] 4
## [2868,] 8
## [2869,] 7
## [2870,] 3
## [2871,] 5
## [2872,] 5
## [2873,] 7
## [2874,] 6
## [2875,] 2
## [2876,] 2
## [2877,] 3
## [2878,] 3
## [2879,] 8
## [2880,] 8
## [2881,] 6
## [2882,] 4
## [2883,] 7
## [2884,] 8
## [2885,] 7
## [2886,] 7
## [2887,] 6
## [2888,] 4
## [2889,] 6
## [2890,] 6
## [2891,] 8
## [2892,] 3
## [2893,] 3
## [2894,] 8
## [2895,] 6
## [2896,] 6
## [2897,] 5
## [2898,] 3
## [2899,] 5
## [2900,] 4
## [2901,] 8
## [2902,] 7
## [2903,] 5
## [2904,] 6
## [2905,] 9
## [2906,] 7
## [2907,] 7
## [2908,] 5
## [2909,] 5
## [2910,] 3
## [2911,] 3
## [2912,] 0
## [2913,] 5
## [2914,] 1
## [2915,] 3
## [2916,] 8
## [2917,] 8
## [2918,] 5
## [2919,] 6
## [2920,] 5
## [2921,] 7
## [2922,] 7
## [2923,] 7
## [2924,] 8
## [2925,] 6
## [2926,] 5
## [2927,] 3
## [2928,] 5
## [2929,] 8
## [2930,] 4
## [2931,] 6
## [2932,] 7
## [2933,] 4
## [2934,] 5
## [2935,] 6
## [2936,] 3
## [2937,] 5
## [2938,] 3
## [2939,] 6
## [2940,] 7
## [2941,] 8
## [2942,] 6
## [2943,] 4
## [2944,] 5
## [2945,] 6
## [2946,] 5
## [2947,] 3
## [2948,] 8
## [2949,] 6
## [2950,] 6
## [2951,] 6
## [2952,] 5
## [2953,] 7
## [2954,] 7
## [2955,] 8
## [2956,] 6
## [2957,] 7
## [2958,] 3
## [2959,] 5
## [2960,] 3
## [2961,] 6
## [2962,] 5
## [2963,] 3
## [2964,] 7
## [2965,] 8
## [2966,] 7
## [2967,] 9
## [2968,] 9
## [2969,] 8
## [2970,] 8
## [2971,] 5
## [2972,] 4
## [2973,] 2
## [2974,] 7
## [2975,] 7
## [2976,] 8
## [2977,] 4
## [2978,] 7
## [2979,] 7
## [2980,] 1
## [2981,] 7
## [2982,] 2
## [2983,] 4
## [2984,] 2
## [2985,] 8
## [2986,] 5
## [2987,] 6
## [2988,] 9
## [2989,] 4
## [2990,] 7
## [2991,] 5
## [2992,] 5
## [2993,] 2
## [2994,] 6
## [2995,] 4
## [2996,] 1
## [2997,] 3
## [2998,] 6
## [2999,] 6
## [3000,] 4
## [3001,] 4
## [3002,] 5
## [3003,] 8
## [3004,] 7
## [3005,] 8
## [3006,] 3
## [3007,] 4
## [3008,] 4
## [3009,] 5
## [3010,] 6
## [3011,] 3
## [3012,] 3
## [3013,] 9
## [3014,] 6
## [3015,] 7
## [3016,] 7
## [3017,] 5
## [3018,] 7
## [3019,] 5
## [3020,] 6
## [3021,] 3
## [3022,] 7
## [3023,] 5
## [3024,] 3
## [3025,] 4
## [3026,] 4
## [3027,] 3
## [3028,] 7
## [3029,] 7
## [3030,] 7
## [3031,] 8
## [3032,] 7
## [3033,] 8
## [3034,] 5
## [3035,] 2
## [3036,] 3
## [3037,] 6
## [3038,] 8
## [3039,] 6
## [3040,] 4
## [3041,] 7
## [3042,] 3
## [3043,] 5
## [3044,] 6
## [3045,] 5
## [3046,] 7
## [3047,] 6
## [3048,] 8
## [3049,] 6
## [3050,] 7
## [3051,] 6
## [3052,] 5
## [3053,] 5
## [3054,] 4
## [3055,] 4
## [3056,] 7
## [3057,] 6
## [3058,] 7
## [3059,] 1
## [3060,] 1
## [3061,] 6
## [3062,] 6
## [3063,] 6
## [3064,] 3
## [3065,] 6
## [3066,] 6
## [3067,] 7
## [3068,] 7
## [3069,] 3
## [3070,] 6
## [3071,] 7
## [3072,] 8
## [3073,] 7
## [3074,] 7
## [3075,] 9
## [3076,] 7
## [3077,] 6
## [3078,] 6
## [3079,] 9
## [3080,] 7
## [3081,] 7
## [3082,] 4
## [3083,] 4
## [3084,] 5
## [3085,] 2
## [3086,] 7
## [3087,] 4
## [3088,] 6
## [3089,] 8
## [3090,] 9
## [3091,] 6
## [3092,] 5
## [3093,] 7
## [3094,] 4
## [3095,] 5
## [3096,] 6
## [3097,] 4
## [3098,] 8
## [3099,] 7
## [3100,] 3
## [3101,] 3
## [3102,] 1
## [3103,] 9
## [3104,] 7
## [3105,] 7
## [3106,] 4
## [3107,] 7
## [3108,] 4
## [3109,] 7
## [3110,] 5
## [3111,] 5
## [3112,] 7
## [3113,] 6
## [3114,] 9
## [3115,] 5
## [3116,] 5
## [3117,] 5
## [3118,] 6
## [3119,] 7
## [3120,] 8
## [3121,] 4
## [3122,] 9
## [3123,] 6
## [3124,] 3
## [3125,] 2
## [3126,] 6
## [3127,] 6
## [3128,] 7
## [3129,] 4
## [3130,] 4
## [3131,] 5
## [3132,] 3
## [3133,] 6
## [3134,] 8
## [3135,] 4
## [3136,] 7
## [3137,] 5
## [3138,] 5
## [3139,] 5
## [3140,] 6
## [3141,] 5
## [3142,] 7
## [3143,] 4
## [3144,] 1
## [3145,] 4
## [3146,] 4
## [3147,] 5
## [3148,] 5
## [3149,] 8
## [3150,] 8
## [3151,] 5
## [3152,] 5
## [3153,] 5
## [3154,] 3
## [3155,] 6
## [3156,] 5
## [3157,] 4
## [3158,] 3
## [3159,] 4
## [3160,] 3
## [3161,] 8
## [3162,] 4
## [3163,] 9
## [3164,] 3
## [3165,] 6
## [3166,] 9
## [3167,] 6
## [3168,] 5
## [3169,] 6
## [3170,] 8
## [3171,] 7
## [3172,] 5
## [3173,] 8
## [3174,] 8
## [3175,] 4
## [3176,] 7
## [3177,] 6
## [3178,] 2
## [3179,] 2
## [3180,] 6
## [3181,] 6
## [3182,] 8
## [3183,] 7
## [3184,] 6
## [3185,] 7
## [3186,] 5
## [3187,] 8
## [3188,] 7
## [3189,] 7
## [3190,] 7
## [3191,] 4
## [3192,] 4
## [3193,] 5
## [3194,] 3
## [3195,] 4
## [3196,] 7
## [3197,] 6
## [3198,] 5
## [3199,] 6
## [3200,] 8
## [3201,] 4
## [3202,] 7
## [3203,] 9
## [3204,] 7
## [3205,] 6
## [3206,] 7
## [3207,] 5
## [3208,] 7
## [3209,] 5
## [3210,] 2
## [3211,] 6
## [3212,] 3
## [3213,] 6
## [3214,] 5
## [3215,] 7
## [3216,] 6
## [3217,] 6
## [3218,] 5
## [3219,] 7
## [3220,] 4
## [3221,] 3
## [3222,] 6
## [3223,] 5
## [3224,] 2
## [3225,] 5
## [3226,] 7
## [3227,] 6
## [3228,] 2
## [3229,] 7
## [3230,] 4
## [3231,] 3
## [3232,] 2
## [3233,] 3
## [3234,] 8
## [3235,] 6
## [3236,] 9
## [3237,] 3
## [3238,] 8
## [3239,] 6
## [3240,] 4
## [3241,] 4
## [3242,] 6
## [3243,] 6
## [3244,] 7
## [3245,] 8
## [3246,] 2
## [3247,] 9
## [3248,] 7
## [3249,] 6
## [3250,] 6
## [3251,] 4
## [3252,] 8
## [3253,] 8
## [3254,] 7
## [3255,] 4
## [3256,] 7
## [3257,] 6
## [3258,] 7
## [3259,] 5
## [3260,] 8
## [3261,] 6
## [3262,] 6
## [3263,] 6
## [3264,] 9
## [3265,] 6
## [3266,] 6
## [3267,] 7
## [3268,] 7
## [3269,] 5
## [3270,] 4
## [3271,] 4
## [3272,] 5
## [3273,] 5
## [3274,] 6
## [3275,] 9
## [3276,] 8
## [3277,] 7
## [3278,] 3
## [3279,] 1
## [3280,] 9
## [3281,] 8
## [3282,] 8
## [3283,] 8
## [3284,] 7
## [3285,] 7
## [3286,] 6
## [3287,] 2
## [3288,] 5
## [3289,] 5
## [3290,] 5
## [3291,] 7
## [3292,] 7
## [3293,] 8
## [3294,] 4
## [3295,] 5
## [3296,] 5
## [3297,] 5
## [3298,] 6
## [3299,] 6
## [3300,] 6
## [3301,] 7
## [3302,] 8
## [3303,] 6
## [3304,] 8
## [3305,] 5
## [3306,] 4
## [3307,] 4
## [3308,] 7
## [3309,] 7
## [3310,] 4
## [3311,] 4
## [3312,] 5
## [3313,] 4
## [3314,] 4
## [3315,] 4
## [3316,] 4
## [3317,] 7
## [3318,] 6
## [3319,] 5
## [3320,] 7
## [3321,] 6
## [3322,] 3
## [3323,] 8
## [3324,] 4
## [3325,] 7
## [3326,] 5
## [3327,] 5
## [3328,] 5
## [3329,] 5
## [3330,] 5
## [3331,] 5
## [3332,] 7
## [3333,] 6
## [3334,] 4
## [3335,] 8
## [3336,] 7
## [3337,] 7
## [3338,] 6
## [3339,] 7
## [3340,] 6
## [3341,] 8
## [3342,] 5
## [3343,] 6
## [3344,] 7
## [3345,] 5
## [3346,] 4
## [3347,] 1
## [3348,] 8
## [3349,] 8
## [3350,] 5
## [3351,] 4
## [3352,] 5
## [3353,] 5
## [3354,] 3
## [3355,] 1
## [3356,] 6
## [3357,] 5
## [3358,] 6
## [3359,] 6
## [3360,] 5
## [3361,] 9
## [3362,] 4
## [3363,] 5
## [3364,] 7
## [3365,] 4
## [3366,] 6
## [3367,] 6
## [3368,] 7
## [3369,] 6
## [3370,] 6
## [3371,] 7
## [3372,] 9
## [3373,] 9
## [3374,] 8
## [3375,] 8
## [3376,] 9
## [3377,] 4
## [3378,] 6
## [3379,] 5
## [3380,] 6
## [3381,] 8
## [3382,] 5
## [3383,] 4
## [3384,] 7
## [3385,] 5
## [3386,] 3
## [3387,] 5
## [3388,] 4
## [3389,] 5
## [3390,] 3
## [3391,] 8
## [3392,] 7
## [3393,] 3
## [3394,] 8
## [3395,] 8
## [3396,] 8
## [3397,] 7
## [3398,] 9
## [3399,] 8
## [3400,] 8
## [3401,] 5
## [3402,] 6
## [3403,] 8
## [3404,] 8
## [3405,] 7
## [3406,] 9
## [3407,] 5
## [3408,] 5
## [3409,] 7
## [3410,] 5
## [3411,] 7
## [3412,] 5
## [3413,] 6
## [3414,] 4
## [3415,] 2
## [3416,] 3
## [3417,] 3
## [3418,] 9
## [3419,] 4
## [3420,] 7
## [3421,] 8
## [3422,] 6
## [3423,] 7
## [3424,] 7
## [3425,] 7
## [3426,] 8
## [3427,] 7
## [3428,] 9
## [3429,] 7
## [3430,] 9
## [3431,] 9
## [3432,] 7
## [3433,] 9
## [3434,] 7
## [3435,] 8
## [3436,] 8
## [3437,] 4
## [3438,] 5
## [3439,] 5
## [3440,] 4
## [3441,] 6
## [3442,] 5
## [3443,] 4
## [3444,] 7
## [3445,] 6
## [3446,] 7
## [3447,] 7
## [3448,] 1
## [3449,] 5
## [3450,] 5
## [3451,] 6
## [3452,] 7
## [3453,] 7
## [3454,] 7
## [3455,] 8
## [3456,] 8
## [3457,] 8
## [3458,] 7
## [3459,] 7
## [3460,] 5
## [3461,] 9
## [3462,] 7
## [3463,] 9
## [3464,] 1
## [3465,] 2
## [3466,] 3
## [3467,] 4
## [3468,] 2
## [3469,] 5
## [3470,] 3
## [3471,] 3
## [3472,] 3
## [3473,] 5
## [3474,] 5
## [3475,] 7
## [3476,] 8
## [3477,] 8
## [3478,] 8
## [3479,] 3
## [3480,] 5
## [3481,] 5
## [3482,] 6
## [3483,] 8
## [3484,] 5
## [3485,] 3
## [3486,] 8
## [3487,] 6
## [3488,] 7
## [3489,] 7
## [3490,] 6
## [3491,] 9
## [3492,] 2
## [3493,] 5
## [3494,] 4
## [3495,] 6
## [3496,] 5
## [3497,] 5
## [3498,] 3
## [3499,] 3
## [3500,] 7
## [3501,] 7
## [3502,] 7
## [3503,] 7
## [3504,] 6
## [3505,] 9
## [3506,] 6
## [3507,] 6
## [3508,] 4
## [3509,] 8
## [3510,] 5
## [3511,] 9
## [3512,] 5
## [3513,] 8
## [3514,] 7
## [3515,] 4
## [3516,] 7
## [3517,] 5
## [3518,] 4
## [3519,] 4
## [3520,] 4
## [3521,] 7
## [3522,] 5
## [3523,] 7
## [3524,] 8
## [3525,] 8
## [3526,] 7
## [3527,] 2
## [3528,] 5
## [3529,] 4
## [3530,] 6
## [3531,] 5
## [3532,] 2
## [3533,] 4
## [3534,] 7
## [3535,] 9
## [3536,] 6
## [3537,] 8
## [3538,] 6
## [3539,] 6
## [3540,] 8
## [3541,] 5
## [3542,] 8
## [3543,] 3
## [3544,] 5
## [3545,] 6
## [3546,] 2
## [3547,] 5
## [3548,] 8
## [3549,] 7
## [3550,] 8
## [3551,] 9
## [3552,] 8
## [3553,] 7
## [3554,] 7
## [3555,] 7
## [3556,] 5
## [3557,] 7
## [3558,] 6
## [3559,] 8
## [3560,] 6
## [3561,] 8
## [3562,] 8
## [3563,] 9
## [3564,] 7
## [3565,] 8
## [3566,] 5
## [3567,] 5
## [3568,] 5
## [3569,] 4
## [3570,] 7
## [3571,] 6
## [3572,] 7
## [3573,] 4
## [3574,] 8
## [3575,] 8
## [3576,] 7
## [3577,] 7
## [3578,] 3
## [3579,] 7
## [3580,] 8
## [3581,] 8
## [3582,] 5
## [3583,] 8
## [3584,] 8
## [3585,] 4
## [3586,] 8
## [3587,] 5
## [3588,] 6
## [3589,] 6
## [3590,] 3
## [3591,] 5
## [3592,] 4
## [3593,] 7
## [3594,] 5
## [3595,] 2
## [3596,] 3
## [3597,] 7
## [3598,] 8
## [3599,] 7
## [3600,] 3
## [3601,] 9
## [3602,] 4
## [3603,] 6
## [3604,] 6
## [3605,] 8
## [3606,] 8
## [3607,] 2
## [3608,] 4
## [3609,] 6
## [3610,] 6
## [3611,] 2
## [3612,] 8
## [3613,] 4
## [3614,] 5
## [3615,] 5
## [3616,] 1
## [3617,] 7
## [3618,] 5
## [3619,] 4
## [3620,] 6
## [3621,] 7
## [3622,] 7
## [3623,] 5
## [3624,] 4
## [3625,] 8
## [3626,] 5
## [3627,] 7
## [3628,] 9
## [3629,] 8
## [3630,] 8
## [3631,] 3
## [3632,] 6
## [3633,] 7
## [3634,] 6
## [3635,] 4
## [3636,] 5
## [3637,] 5
## [3638,] 6
## [3639,] 4
## [3640,] 6
## [3641,] 4
## [3642,] 7
## [3643,] 5
## [3644,] 5
## [3645,] 6
## [3646,] 6
## [3647,] 6
## [3648,] 7
## [3649,] 7
## [3650,] 7
## [3651,] 7
## [3652,] 9
## [3653,] 7
## [3654,] 8
## [3655,] 8
## [3656,] 8
## [3657,] 8
## [3658,] 8
## [3659,] 5
## [3660,] 5
## [3661,] 6
## [3662,] 3
## [3663,] 8
## [3664,] 3
## [3665,] 4
## [3666,] 6
## [3667,] 5
## [3668,] 7
## [3669,] 7
## [3670,] 8
## [3671,] 5
## [3672,] 4
## [3673,] 4
## [3674,] 5
## [3675,] 8
## [3676,] 3
## [3677,] 6
## [3678,] 9
## [3679,] 5
## [3680,] 1
## [3681,] 9
## [3682,] 7
## [3683,] 6
## [3684,] 4
## [3685,] 7
## [3686,] 6
## [3687,] 8
## [3688,] 5
## [3689,] 6
## [3690,] 5
## [3691,] 9
## [3692,] 5
## [3693,] 6
## [3694,] 7
## [3695,] 7
## [3696,] 8
## [3697,] 6
## [3698,] 4
## [3699,] 6
## [3700,] 7
## [3701,] 6
## [3702,] 4
## [3703,] 5
## [3704,] 5
## [3705,] 5
## [3706,] 6
## [3707,] 7
## [3708,] 3
## [3709,] 6
## [3710,] 8
## [3711,] 9
## [3712,] 9
## [3713,] 0
## [3714,] 2
## [3715,] 8
## [3716,] 6
## [3717,] 7
## [3718,] 9
## [3719,] 6
## [3720,] 7
## [3721,] 9
## [3722,] 5
## [3723,] 7
## [3724,] 5
## [3725,] 8
## [3726,] 4
## [3727,] 5
## [3728,] 4
## [3729,] 8
## [3730,] 9
## [3731,] 7
## [3732,] 6
## [3733,] 8
## [3734,] 9
## [3735,] 4
## [3736,] 5
## [3737,] 6
## [3738,] 7
## [3739,] 6
## [3740,] 4
## [3741,] 5
## [3742,] 6
## [3743,] 6
## [3744,] 3
## [3745,] 2
## [3746,] 6
## [3747,] 5
## [3748,] 7
## [3749,] 6
## [3750,] 9
## [3751,] 9
## [3752,] 7
## [3753,] 6
## [3754,] 7
## [3755,] 8
## [3756,] 7
## [3757,] 9
## [3758,] 7
## [3759,] 7
## [3760,] 7
## [3761,] 6
## [3762,] 6
## [3763,] 8
## [3764,] 5
## [3765,] 7
## [3766,] 5
## [3767,] 7
## [3768,] 6
## [3769,] 4
## [3770,] 4
## [3771,] 8
## [3772,] 7
## [3773,] 4
## [3774,] 5
## [3775,] 7
## [3776,] 9
## [3777,] 6
## [3778,] 6
## [3779,] 9
## [3780,] 9
## [3781,] 6
## [3782,] 7
## [3783,] 3
## [3784,] 4
## [3785,] 2
## [3786,] 2
## [3787,] 4
## [3788,] 4
## [3789,] 4
## [3790,] 5
## [3791,] 3
## [3792,] 5
## [3793,] 6
## [3794,] 6
## [3795,] 6
## [3796,] 5
## [3797,] 5
## [3798,] 6
## [3799,] 5
## [3800,] 4
## [3801,] 8
## [3802,] 8
## [3803,] 8
## [3804,] 7
## [3805,] 6
## [3806,] 3
## [3807,] 6
## [3808,] 2
## [3809,] 7
## [3810,] 3
## [3811,] 3
## [3812,] 4
## [3813,] 6
## [3814,] 7
## [3815,] 6
## [3816,] 3
## [3817,] 6
## [3818,] 5
## [3819,] 5
## [3820,] 3
## [3821,] 5
## [3822,] 8
## [3823,] 0
## [3824,] 7
## [3825,] 7
## [3826,] 8
## [3827,] 4
## [3828,] 6
## [3829,] 5
## [3830,] 3
## [3831,] 3
## [3832,] 3
## [3833,] 2
## [3834,] 7
## [3835,] 7
## [3836,] 7
## [3837,] 5
## [3838,] 8
## [3839,] 6
## [3840,] 5
## [3841,] 3
## [3842,] 7
## [3843,] 6
## [3844,] 4
## [3845,] 6
## [3846,] 8
## [3847,] 6
## [3848,] 7
## [3849,] 7
## [3850,] 5
## [3851,] 5
## [3852,] 2
## [3853,] 5
## [3854,] 5
## [3855,] 6
## [3856,] 8
## [3857,] 8
## [3858,] 8
## [3859,] 9
## [3860,] 8
## [3861,] 7
## [3862,] 9
## [3863,] 4
## [3864,] 7
## [3865,] 7
## [3866,] 7
## [3867,] 6
## [3868,] 8
## [3869,] 8
## [3870,] 6
## [3871,] 7
## [3872,] 8
## [3873,] 6
## [3874,] 8
## [3875,] 5
## [3876,] 5
## [3877,] 7
## [3878,] 9
## [3879,] 6
## [3880,] 5
## [3881,] 6
## [3882,] 6
## [3883,] 4
## [3884,] 2
## [3885,] 6
## [3886,] 9
## [3887,] 6
## [3888,] 8
## [3889,] 8
## [3890,] 6
## [3891,] 7
## [3892,] 8
## [3893,] 6
## [3894,] 6
## [3895,] 7
## [3896,] 5
## [3897,] 3
## [3898,] 7
## [3899,] 7
## [3900,] 7
## [3901,] 9
## [3902,] 6
## [3903,] 6
## [3904,] 7
## [3905,] 7
## [3906,] 6
## [3907,] 2
## [3908,] 7
## [3909,] 5
## [3910,] 3
## [3911,] 7
## [3912,] 6
## [3913,] 7
## [3914,] 6
## [3915,] 7
## [3916,] 7
## [3917,] 6
## [3918,] 8
## [3919,] 7
## [3920,] 5
## [3921,] 7
## [3922,] 7
## [3923,] 5
## [3924,] 7
## [3925,] 8
## [3926,] 6
## [3927,] 4
## [3928,] 6
## [3929,] 7
## [3930,] 4
## [3931,] 4
## [3932,] 3
## [3933,] 7
## [3934,] 6
## [3935,] 8
## [3936,] 5
## [3937,] 4
## [3938,] 6
## [3939,] 7
## [3940,] 5
## [3941,] 6
## [3942,] 3
## [3943,] 5
## [3944,] 3
## [3945,] 5
## [3946,] 7
## [3947,] 4
## [3948,] 8
## [3949,] 2
## [3950,] 4
## [3951,] 7
## [3952,] 8
## [3953,] 7
## [3954,] 5
## [3955,] 3
## [3956,] 6
## [3957,] 7
## [3958,] 4
## [3959,] 3
## [3960,] 8
## [3961,] 5
## [3962,] 7
## [3963,] 7
## [3964,] 8
## [3965,] 6
## [3966,] 5
## [3967,] 3
## [3968,] 9
## [3969,] 5
## [3970,] 2
## [3971,] 2
## [3972,] 4
## [3973,] 1
## [3974,] 4
## [3975,] 5
## [3976,] 8
## [3977,] 4
## [3978,] 4
## [3979,] 5
## [3980,] 9
## [3981,] 4
## [3982,] 5
## [3983,] 5
## [3984,] 7
## [3985,] 4
## [3986,] 9
## [3987,] 5
## [3988,] 7
## [3989,] 7
## [3990,] 4
## [3991,] 7
## [3992,] 0
## [3993,] 4
## [3994,] 5
## [3995,] 5
## [3996,] 7
## [3997,] 6
## [3998,] 9
## [3999,] 8
## [4000,] 7data.frame(obs = ppd) %>%
ggplot(aes(x=obs)) +
geom_histogram() +
labs(x="Observed water (out of 9)")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Sampling parameters from the prior
Oops, we’ve jumped ahead of ourselves! Best practice is to simulate values from your prior first and check to see if those priors are reasonable.
m1p <-
brm(data = list(w = 6),
family = binomial(link = "identity"),
w | trials(9) ~ 0 + Intercept,
prior(beta(1, 1), class = b, lb = 0, ub = 1),
iter = 5000, warmup = 1000, seed = 3, chains=1,
sample_prior = "only")samples_from_prior = as_draws_df(m1p)
samples_from_prior %>%
ggplot(aes(x=b_Intercept)) +
geom_density(fill = "grey", color = "white") +
labs(x="Proportion water", title="Prior")
Simulating observations from the prior
We may also want the PRIOR PREDICTIVE DISTRIBUTION which is the expected observiations given our prior.
prior_pd = posterior_predict(m1p)
data.frame(obs = prior_pd) %>%
ggplot(aes(x=obs)) +
geom_histogram() +
labs(x="Observed water (out of 9)")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Simulating from your priors – prior predictive simulation – is an essential part of modeling. This allows you to see what your choices imply about the data. You’ll be able to diagnose bad choices.
an aside about learning in R
At this point in the course, I’m going to start throwing a lot of code at you. Do I expect you to memorize this code? Of course not.
Do you need to understand every single thing that’s happening in the code? Nope.
But, you’ll learn a lot by taking the time to figure out what’s happening in a code chunk. Class time will frequently include exercises where I ask you to adapt code I’ve shared in the slides to a new dataset or to answer a new problem. When doing so, go back through the old code and figure out what’s going on. Run the code one line at a time. Always observe the output and take some time to look at the object that was created or modified. Here are some functions that will be extremely useful:
str() # what kind of object is this? what is its structure?
dim() # what are the dimensions (rows/columns) of this object
head() # give me the first bit of this objectstr(prior_pd)## int [1:4000, 1] 6 0 1 1 9 7 0 5 4 5 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : NULLdim(prior_pd)## [1] 4000 1head(prior_pd)## [,1]
## [1,] 6
## [2,] 0
## [3,] 1
## [4,] 1
## [5,] 9
## [6,] 7Continous outcomes
The globe tossing example is cute and easy to work with, but let’s move towards the kinds of variables we more frequently work with. Let’s create a model for some outcome, \(y\) that is continuous.
\[\begin{align*} y_i &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(0, 10) \\ \sigma &\sim \text{Uniform}(0, 5) \end{align*}\]
set.seed(9)
y = rnorm(n = 31, mean = 4, sd = .5)
m2 = brm(
data = list(y=y),
family = gaussian,
y ~ 1,
prior = c(prior( normal(0,10), class=Intercept),
prior( uniform(0,5), class=sigma)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
file = here("files/models/m21.2")
)An example: weight and height
Using the Howell data (don’t load the rethinking package
because it can interfere with brms).
data("Howell1", package = "rethinking")
d <- Howell1
str(d)## 'data.frame': 544 obs. of 4 variables:
## $ height: num 152 140 137 157 145 ...
## $ weight: num 47.8 36.5 31.9 53 41.3 ...
## $ age : num 63 63 65 41 51 35 32 27 19 54 ...
## $ male : int 1 0 0 1 0 1 0 1 0 1 ...library(measurements)
d$height <- conv_unit(d$height, from = "cm", to = "feet")
d$weight <- conv_unit(d$weight, from = "kg", to = "lbs")
rethinking::precis(d)## mean sd 5.5% 94.5% histogram
## height 4.5362072 0.9055921 2.661042 5.4375 ▁▁▁▁▂▂▇▇▁
## weight 78.5079631 32.4502290 20.636856 120.1583 ▁▂▃▂▂▁▁▃▅▇▇▃▂▁
## age 29.3443934 20.7468882 1.000000 66.1350 ▇▅▅▃▅▂▂▁▁
## male 0.4724265 0.4996986 0.000000 1.0000 ▇▁▁▁▁▁▁▁▁▇d2 <- d[ d$age >= 18, ]exercise
Write a mathematical model for the weights in this data set. (Don’t worry about predicting from other variables yet.)
solution
\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
\[\begin{align*} w &\sim \text{Normal}(\mu, \sigma) \\ \mu &\sim \text{Normal}(130, 20) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
exercise
Simulate from your priors (parameters values and prior predictive values).
solution
Sample from your priors:
m3p = brm(
data = d2,
family = gaussian,
weight ~ 1,
prior = c(prior( normal(130,20), class=Intercept),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
sample_prior = "only")Sampling parameter estimates for the Intercept.
pairs(m3p)## Warning: Only one chain in 'x'. This plot is more useful with multiple chains.
This is a different (shorter) way to plot your posterior. Good things: it automatically includes the scatterplot so you can see the implications of how these parameters correlate. Bad things: not customizable and not useable when you have a lot of parameters.
Simulate values of weight.
prior_pd = posterior_predict(m3p)
dim(prior_pd)## [1] 4000 352as.data.frame(prior_pd) %>%
pivot_longer(everything()) %>%
ggplot(aes(x=value)) +
geom_histogram() +
labs(x="Expected observed weights (based on prior)")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Another shorter way:
pp_check(m3p)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
Fit the model
m3 = brm(
data = d2,
family = gaussian,
weight ~ 1,
prior = c(prior( normal(130,20), class=Intercept),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
file = here("files/models/m21.3"))posterior_summary(m3)## Estimate Est.Error Q2.5 Q97.5
## b_Intercept 99.221909 0.75720356 97.751074 100.737488
## sigma 14.291987 0.55198322 13.254363 15.428183
## Intercept 99.221909 0.75720356 97.751074 100.737488
## lprior -8.318377 0.05827127 -8.433538 -8.203915
## lp__ -1441.294276 1.02695551 -1444.235605 -1440.292326pairs(m3)## Warning: Only one chain in 'x'. This plot is more useful with multiple chains.
posterior_predict(m3) %>%
as.data.frame() %>%
pivot_longer(everything()) %>%
ggplot(aes(x=value)) +
geom_density(fill = "grey", color = "white") +
geom_density( aes(x = weight), data=d2, inherit.aes = F) 
pp_check(m3)## Using 10 posterior draws for ppc type 'dens_overlay' by default.
Adding in a linear component
We might assume that height and weight are associated with each other. Indeed, within our sample:
plot(d2$weight ~ d2$height)
exercise
Update your mathematical model to incorporate height.
\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta h_i \\ \alpha &\sim \text{Normal}(??, ??) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
\(=\) is deterministic – once we know other variables, \(\mu_i\) is known with certainty
made-up parameters are the targets of learning
exercise
Update your mathematical model to incorporate height.
\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
\(=\) is deterministic – once we know other variables, \(\mu_i\) is known with certainty
made-up parameters are the targets of learning
\[\begin{align*} w_i &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta (h_i - \bar{h}) \\ \alpha &\sim \text{Normal}(130, 20) \\ \beta &\sim \text{Normal}(0, 25) \\ \sigma &\sim \text{Uniform}(0, 25) \\ \end{align*}\]
To update our brms code:
d2$height_c = d2$height - mean(d2$height)
m4p = brm(
data = d2,
family = gaussian,
weight ~ 1 + height_c,
prior = c(prior( normal(130,20), class=Intercept),
prior( normal(0,25), class=b),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
sample_prior = "only")set.seed(9)
samples_from_prior = as_draws_df(m4p)
str(samples_from_prior)## draws_df [4,000 × 9] (S3: draws_df/draws/tbl_df/tbl/data.frame)
## $ b_Intercept: num [1:4000] 167 155.6 149.6 88.7 176.4 ...
## $ b_height_c : num [1:4000] -57.9 -57.2 -18.5 52 -47.3 ...
## $ sigma : num [1:4000] 14.1 11.8 13.3 2.1 21.1 ...
## $ Intercept : num [1:4000] 167 155.6 149.6 88.7 176.4 ...
## $ lprior : num [1:4000] -15.7 -14.7 -12 -15.6 -15.8 ...
## $ lp__ : num [1:4000] -13.8 -12.9 -10.2 -14.9 -14.6 ...
## $ .chain : int [1:4000] 1 1 1 1 1 1 1 1 1 1 ...
## $ .iteration : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...
## $ .draw : int [1:4000] 1 2 3 4 5 6 7 8 9 10 ...d2 %>%
ggplot(aes(x=height_c, y=weight)) +
geom_blank() +
geom_abline( aes(intercept=b_Intercept, slope=b_height_c),
data=samples_from_prior[1:50, ],
alpha=.3) +
scale_x_continuous(name = "height(feet)",
breaks=seq(4,6,by=.5)-mean(d2$height),
labels=seq(4,6,by=.5))
Describe in words what’s wrong with our priors.
Slope should not be negative. How can we fix this?
Could use a uniform distribution bounded by 0.
m4p = brm(
data = d2,
family = gaussian,
weight ~ 1 + height_c,
prior = c(prior( normal(130,20), class=Intercept),
prior( uniform(0,25), class=b),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
sample_prior = "only")
exercise
Fit the new weight model to the data.
solution
m4 = brm(
data = d2,
family = gaussian,
weight ~ 1 + height_c,
prior = c(prior( normal(130,20), class=Intercept),
prior( uniform(0,25), class=b),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
file=here("files/models/m21.4"))
posterior_summary(m4)## Estimate Est.Error Q2.5 Q97.5
## b_Intercept 99.21508 0.5154139 98.17317 100.20955
## b_height_c 24.74661 0.2609169 24.05150 24.99355
## sigma 10.36090 0.3898531 9.65234 11.12228
## Intercept 99.21508 0.5154139 98.17317 100.20955
## lprior -11.53739 0.0396881 -11.61861 -11.46176
## lp__ -1332.15882 1.3936194 -1335.90639 -1330.52002exercise
Draw lines from the posterior distribution and plot with the data.
solution
set.seed(9)
samples_from_posterior = as_draws_df(m4)
d2%>%
ggplot(aes(x=height_c, y=weight)) +
geom_point(size=.5) +
geom_abline( aes(intercept=b_Intercept, slope=b_height_c),
data=samples_from_posterior[1:50, ],
alpha=.3,
color="#1c5253") +
scale_x_continuous(name = "height(feet)",
breaks=seq(4,6,by=.5)-mean(d2$height),
labels=seq(4,6,by=.5))
A side note: a major concern or critique of Bayesian analysis is that the subjectivity of the priors allow for nefarious behavior. “Putting our thumbs on the scale,” so to speak. But priors are quickly overwhelmed by data. Case in point:
m4e = brm(
data = d2,
family = gaussian,
weight ~ 1 + height_c,
prior = c(prior( normal(130,20), class=Intercept),
prior( normal(-5,5), class=b),
prior( uniform(0,25), class=sigma, lb=0, ub=25)),
iter = 5000, warmup = 1000, seed = 3, chains=1,
file=here("files/models/m21.4e"))
posterior_summary(m4e)## Estimate Est.Error Q2.5 Q97.5
## b_Intercept 99.197733 0.5035653 98.214092 100.15932
## b_height_c 35.775818 1.8832462 32.177572 39.50070
## sigma 9.529429 0.3687178 8.839379 10.27233
## Intercept 99.197733 0.5035653 98.214092 100.15932
## lprior -44.172476 3.0738968 -50.456409 -38.49994
## lp__ -1334.629214 1.1849503 -1337.645046 -1333.23509You’ll only really get into trouble with uniform priors that have a boundary, if true population parameter is outside your boundary. A good rule of thumb is to avoid the uniform distribution. We’ll cover other options for priors for \(\sigma\) in future lectures, but as a preview, the exponential distribution works very well for this!